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Abstract Our understanding of the elasticity and rhe- 
ology of disordered materials, such as granular piles, 
foams, emulsions or dense suspensions relies on improv- 
ing experimental tools to characterize their behaviour 
at the particle scale. While 2D observations are now 
routinely carried out in laboratories, 3D measurements 
remain a challenge. In this paper, we use a simple model 
system, a packing of soft elastic spheres, to illustrate the 
capability of X-ray microtomography to characterise 
the internal structure and local behaviour of granular 
systems. Image analysis techniques can resolve grain 
positions, shapes and contact areas; this is used to in- 
vestigate the material's microstructure and its evolu- 
tion upon strain. In addition to morphological mea- 
surements, we develop a technique to quantify contact 
forces and estimate the internal stress tensor. As will 
be illustrated in this paper, this opens the door to a 
broad array of static and dynamical measurements in 
3D disordered systems. 
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Despite the practical importance of macroscopic dis- 
ordered materials such as soil, foams and emulsions the 
rules that govern their mechanical behaviour remain 
poorly understood. In the case of granular materials, 
these studies are essential to the advancement of re- 
lated industrial processes and also to the prediction 
of often catastrophic geological phenomena (avalanches 
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of slurries, earthquakes). Although empirical constitu- 
tive equations are nowadays able to partly model the 
response of these systems, a unified multiscale frame- 
work does not exist yet. Over the past 20 years, foams 
and granular systems have emerged as a model system 
for low temperature glasses [T]. Many theories and ex- 
periments have been designed to address the relation- 
ship between the internal structure and macroscopic 
behaviour of such systems. These theories and models 
often use approaches derived from statistical physics [2j 
[3] . Although the behaviour of foams seems today rather 
well understood [4], the complexity of granular mat- 
ters behaviour at the microscopic scale has prevented 
a proper physical model of granular systems. A large 
part of the problem arises from the strongly nonlinear 
contact law between rigid bodies coupled with the dy- 
namical nature of the contact network. In this context, 
the spatial distribution and temporal evolution of the 
force network has become a highly sought after quantity 
that expands from micro scale (grain-grain contacts) to 
macro scale (granular assembly). 

Early simulations have provided important insights 
on the internal force distributions of a granular pile [5j 
EJ[7]. Two-dimensional (2D) experiments managed to 
keep up with the pace of advancement in simulations 
[8,9,10 , 3D experiments however, have been limited to 
measurements of the distribution of forces only at the 
boundaries of the containers jTTHT2l fT3| rT4] . These ex- 
perimental methods could not have access to the spa- 
tial arrangement of the contact force network in the 
bulk of granular assemply. Moreover they were unable 
to determine structural features such as force chains 
and arching which have been postulated as the signa- 
ture of jamming [T5] . 

In recent years, a range of tools have been devel- 
oped to apprehend the 3D nature of bulk properties. 
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Fig. 1 Sketch of the compression cell and consecutive stages of 
compression. 



X-ray or MRI techniques have allowed to study dy- 
namic properties of granular systems such as flow pro- 
files and shear banding, at a mesoscopic scale [T6 | fT7 | 
18,19 . 3D reconstructions of compressed emulsion sys- 
tems using confocal microscopy have provided the first 
measurements of the bulk force distributions [20 ,21] by 
estimating the contact force from contact geometry. A 
similar technique has also been used to characterise spa- 
tial correlations of forces inside 3D piles of frictionless 
liquid droplets [22]. Although these observations are 
very valuable, the systems used are significantly dif- 
ferent from real granular pile which has very different 
contact laws and frictional properties. 



X-ray computed tomography has been used for static 
characterisation of large packings of spherical grains 
[23,24,25 r In this study, we apply this technique to 
a model system made of soft millimetric elastic beads. 
In addition to the measurement of traditional struc- 
tural quantities (packing fraction, coordination num- 
ber etc.), the contact force is measured from the con- 
tact area between grains which can be integrated at a 
mesoscopic scale to estimate the local stress tensor in 
the pile. We then apply a series of external forces to 
our model system and investigate its response to these 
controlled loadings. This approach delivers results con- 
sistent with the existing literature and is promising as 
a generic tool to study the local, non-linear mechanics 
of granular assemblies. 



1 Experiment and methodology 

1.1 Experimental setup 

The packing studied in this paper is made of relatively 
monodisperse spherical rubber balls of diameter 3.10 ± 
0.05 mmB The maximum and minimum ball diameter 
measured in the packing are 3.19mm and 2.85mm re- 
spectively. The grains are made of commercial rubber 
with a shear modulus estimated at 850kPa . 

A cylindrical PMMA container, with internal diam- 
eter of 44mm is used as a compression cell (see figure 
[l}. The inner wall of the cylinder is lubricated with 
canola oil to reduce the friction. 2020 rubber elastic 
balls are then poured into this chamber. The latter is 
then closed at the top by a piston of diameter slightly 
smaller than the container; it can therefore move freely 
without touching the PMMA cylinder and without al- 
lowing beads to leave the container. A horizontal plat- 
form is rigidly connected by a shaft to the piston so 
that a load can be applied to the sample by placing 
weights on it. A pair of springs connecting the piston 
to the container is used to make sure that, in the ab- 
sence of a load, the piston does not fall down due to its 
own weight. 

The compression cell is attached to a motorised ro- 
tation stage located between a high-resolution microfo- 
cus X-ray source (80kV accelerating voltage and 200/iA 
beam current and a CCD detector of the size 67mm x 
67mm, 2048 x 2048 pixels of size 33.6/im each [561127]. 
The compression cell is rotated by about 0.125 degree 
increments around its vertical axis and a radiographic 
projection of the packing is taken after each rotation 
with an exposure time of 18 s. A total of 2880 projec- 
tions are taken for a complete rotation of the specimen. 
It takes approximately 15 hours to scan the volume 
for each case. After the completion of each scan, tomo- 
graphic cone beam reconstructions are performed on 
these 2880 projections using the a Feldkamp algorithm 
[28] . Tomograms of about 2000 3 voxels are obtained at 
27 microns resolutions. 

We have recorded and analysed the internal geome- 
try of the packing under three different loading condi- 
tions. 

Stage 0: Pre-compression stage with no extra loading; 
the piston is held in balance just above the top of 
the packing without making contact with it (see Fig. 
[I]). Rubber balls are at rest, only bearing their own 
weight (41. lg). The combined weight of the piston 
and the platform (600g) is fully balanced by the 

1 Styrene-Butadine-Rubber (SBR); purchases from Mid- 
Atlntic Rubber Co., USA 



3 



0.01 1 r 



(a) 




(b) (c) 



Fig. 2 (a) Normalized X-ray density histogram of the full image 
volume of the rubber ball sample, (b) Gray scale X-ray density 
map of a slice of rubber ball sample, (c) The same slices after 
phase separation. 

pull of the springs so that at this point there is no 
external force being exerted onto the pile. 

Stage 1: The platform is loaded with 150#(1.47iV). By 
taking into account the opposite force applied by 
the springs, a net weight of 115g is applied to the 
packing, resulting in a compressive strain of e = 
3.78%. We define the axial strains as the ratio of 
the displacement of piston to the initial height of 
the packing (engineering strains). 

Stage 2: A total of 900#(8.837V) is placed on the plat- 
form resulting in a net 780g weight onto the packing. 
The total axial strain at that stage is e = 7.83%. 

1.2 3D Image analysis: segmentation and partitioning 

The tomographic image consists of a cubic array of 
reconstructed linear X-ray attenuation coefficient val- 
ues, each corresponding to a finite volume cube (voxel) 
of the sample. The first step in analysing this data is 
to differentiate the attenuation map into distinct pore 
and grain phases. Ideally, one would wish to have a 
bi-modal distribution giving unambiguous phase sepa- 
ration of the pore and solid phase peaks. This simple 
phase extraction is possible in our rubber ball pack. 
The intensity histogram of the tomogram (Fig. [5|a)) 
shows two distinct peaks associated with the two phases 
(beads and air). The peak centred around 9000 is asso- 
ciated with the beads and the lower peak around 4000 
is associated with the pore phase. 

The segmentation algorithm [29] uses both the at- 
tenuation value and its gradient to detect interfaces 




Fig. 3 Different stages of the watershed algorithm, (a) Two over- 
lapping discs, (b) EDT of the two overlapping discs, (c) Identify- 
ing the local maxima (seeds), (d) and (e) Regions grow from the 
initial seeds by having pixels on their boundary added to them, 
(f) Discs separated. 



between the two phases (grain- void). In particular, a 
sharp gradient of attenuation is required to detect an 
interface. This method, after adjustment [30], provides 
a very robust detection of grain boundaries, but is of 
course limited to void spaces larger than the voxel size. 
In particular, if the gap between two grains is very small 
(order of a voxel), the density gradient across the gap 
between these two grains will not appear as steep as a 
typical grain- void interface. This can result in the detec- 
tion of false contact regions in the segmented datasets 
which consequently gives rise to the following undesired 
morphological effects: i) detection of false contacts be- 
tween grains, and ii) the surface area of real contacts 
appear larger than their actual values. As detailled later 
in this paper, this effect can be corrected in spherical 
grain piles by offsetting the apparent contact area of 
the grains. 

The next step is to reconstruct and label individual 
grains from the geometry of the solid phase. The basic 
assumption in the grain identification algorithm is that 
the boundaries between grains which are not isolated 
coincide with the watershed surfaces of the Euclidean 
distance map of grains (distance to the nearest grain 
boundary) [31]. The entire grain space can be thought 
of as the union of spheres centred on every grain voxel. 
Each sphere radius is given by the Euclidean distance 
value of the voxel at its centre. The next stage is iden- 
tifying all the voxels that are not covered by any larger 
Euclidean spheres. 

Each one of these voxels, which are at the maxima 
of the distance function in their local neighbourhood, 
then forms a seed that will grow into a single grain in 
the following stage of the algorithm. Fig. [3] illustrates a 
simple 2D case where 2 overlapping discs are separated. 
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Fig. 4 A subvolume of the elastic ball compaction experiment. 
The contact zone between touching grains (dark patches) pro- 
vides a measurement of contact force. 



The seed regions essentially grow by having voxels on 
their boundary added to them. Voxels that lie on the 
boundaries of the regions are processed in reverse Eu- 
clidean distance order, i.e. voxels with high distance 
values are processed first. When a voxel is processed, 
it is assigned to the region on whose boundary it lies, 
or, if it lies on more than one region boundary, the 
region whose boundary it first became part of. At the 
end of the algorithm, the grain space will be partitioned 
into grains whose boundaries lie on the watershed sur- 
faces of the Euclidean distance function [25]. The wa- 
tershed based grain partitioning algorithm is compu- 
tationally expensive. Hence it is parallelised (using an 
implementation of the "time warp" discrete event sim- 
ulation protocol [32 J, so it can be used to analyse very 
large datasets (~ 2000 3 ). 



2 Measurements 

The 3D watershed algorithm provides us with large 
amounts of information about the pile structure. Each 
grain has been labelled, its location and shape are known, 
its direct neighbours have been identified, and the ge- 
ometry of the contact they share can been extracted as 
well. This data can now be processed to extract physical 
and mechanical properties. In this section, we present 
a few methods implemented in the context of granular 
systems. These can be separated into structural charac- 
teristics and mechanical measurements. The latter in- 
volves kinematic and force measurements that rely on 
prior knowledge of a contact law for individual grains. 
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Fig. 5 Autocorrelation function of the packing at different stages 
of compression normalized by the volume of a single grain and 
the packing fraction. Note the shift in the peaks as the pressure 
increases. 



2.1 Structural quantities 

2.1.1 Packing fraction 

One of the most readily measurable quantities from 3D 
images is the packing fraction. Packing fraction or ap- 
parent density is defined as the ratio of the volume of 
balls to the total volume. Table [l] summarises the pack- 
ing fractions measured directly from the digital images. 
The packing fractions we measure on the initial pile 
is 57%, which is in the lower end of the spectrum for 
monodisperse beads. Such a density is achievable in our 
system due to the large frictional coefficient of rubber- 
rubber contacts. As the loading is increased, the volume 
fraction increases up to 62%. 

2.1.2 Radial Distribution Function 

From the shapes of the individual grains calculated 
by the Watershed algorithm, the coordinates of each 
grain's centre of mass can be accurately measured. This 
data is in particular required to characterize the mi- 
crostructure of granular systems and to calculate the 
autocorrelation function of the grain centres. 

The radial distribution function (RDF) is a measure 
of the degree of separation of grain centres and their 



Compression stage 
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Loading [gr] 





115 


780 


Strain [%] 





3.78 


7.83 


Packing fraction 


0.57 


0.59 


0.62 


Avg. coord, no. 
raw 


6.28 ± 1.64 


6.63 ± 1.58 


7.13 ± 1.49 


Avg. coord, no. 
filtered 


4.99 ± 1.80 


5.32 ± 1.66 


5.92 ± 1.56 



Table 1 Details of the compression progression. 
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Fig. 6 Histogram of coordination number as a function of load- 
ing for different threshold values: (a) No threshold is applied (us- 
ing the raw data after watershed algorithm), (b) After a threshold 
value of 0.35 mm 2 is subtracted from all contact surface areas. As 
a result the average coordination number is reduced in all three 



density at a given distance. It is calculated by counting 
the number of grains N that are separated from a given 
sphere by a distance in the interval [r, r + dr] , where r 
is the distance between grain centres. For large r the 
number density of the grain centres found in the interval 
[r, r+dr] approaches the average density of grains in the 
packing, p, where p is the packing fraction and R is 
average grain radius. In our calculations, we normalize 
RDF by the average density. Figure [5] shows the RDF in 
our packings throughout compression progression. The 
RDF of the system before insertion of any external force 
(stage 0) shows a prominent peak at r w D g where D g 
is the average grain diameter (D g = 2 x R « 3.10mm). 
The second peak appears at r « 1.95D g and a sub- 
peak approximately at r = V3D g . As the vertical load 
increases (stage 1 and stage 2), the peak of the RDF 
widens and shifts to the left due to the compression 
of touching spheres, hence shortening the centre-centre 
distance between them (see Fig. [5| left)). 

The magnitude of the shift is about 1% and 2% 
for stages 1 and 2. This is to be compared with the 
macroscopic axial strain which is of the other of 4% 
and 8% without any deformation along the order di- 




Fig. 7 (a) Network representation the full packing, (b) An in- 
ternal view of the packing illustrating the connecting grains 



rections. This suggests that the internal compression is 
more isotropic than the loading. 

2.1.3 Mean coordination 

The grain partitioning also provides us with the list 
of contact areas between grains, from which it is pos- 
sible to analyze the contact network of grains. Using 
first the raw data from the watershed algorithm, we 
have measured the coordination (number of contacts) 
of each grain and plotted its distribution for each stage 
of compression (see Figure^ a)). This plot shows a sig- 
nificant increase of grain coordination as the loading 
increases. It is notable that as pressure increases, the 
distributions move to higher values and get slightly nar- 
rower (see Table [TJ. This increase in the coordination 
number is achieved by re-organisation of grains in the 
packing when they are compressed and also grain com- 
pression/deformation process which reduces the overall 
grain-grain distance. A relatively large number of grains 
lose their contacts with some of their immediate neigh- 
bours while almost the same amount gain new contacts 
as the compression progresses. We have shown the his- 
togram of lost/gained contacts in the inset of figure 
[6] Negative values on the horizontal axis represent lost 
contacts while positive values show the number of newly 
gained contacts. Nearly 50% of grains retain their co- 
ordination number during the compression without any 
changes. 

the inherent finite reso- 



As discussed in section 1.2 



lution of the CT leads to systematic bias in fine details 
such as grain contacts; if the gap between two grains is 
in the order of, or less than, the voxel size, grains will 
appear to be in contact in the segmented data. To cor- 
rect for this intrinsic resolution limitation, we offset all 
the contact areas by a small amount, which essentially 
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Fig. 8 Displacement field presented as vector fields (velocity 
fields) for stage 1 (a) and 2 (b) of the experiment. 



corresponds to the apparent contact area of two touch- 



ing perfectly spherical grains (see section 2.2.3 for de- 
tailed discussion). All contacts whose apparent surface 
area is below this threshold, are therefore discarded. 
Figure j6^b) shows a re-plot of figure [6^ a) after applying 
the offset. The isostatic limit for mechanically stable 
packing of spheres in 3D suggests a connectivity of 4 
for frictional and 6 for frictionless systems. In our mea- 
surements, the average coordination number is c± 5 for 
stage (see Figj6^b)) which is in agreement with pre- 
vious measurements [33| l24] . 

From the knowledge of the whole contact network, 
a large number of other morphological and topologi- 
cal quantities can be calculated, such as spatial corre- 
lations in contact orientations, or contact anisotropy. 
Figure [7] shows a reconstruction of the contact network 
through which forces propagate. The statistical proper- 
ties of such networks will be developed in further stud- 
ies. In what follows, we will focus on the determination 
of mechanical forces within the pile, at the micro scale, 
as well as mesoscopic scale. 



2.2 Kinematics and dynamics 
2.2.1 Displacement fields 

We calculate next the coordinates of the centroid of the 
grains for all three scans and then track each individual 
grain throughout the experiment. As an illustration, we 
render the 3D vector field derived from tracking grains 
during the compression for both stages of compression, 
see figures |8ja-b). 

To gain a qualitative insight into the displacement 
field, we measure the average displacement of grain cen- 
tres in both longitudinal and radial directions. Figure 
[9^a) demonstrates such displacement along the load- 
ing axis. As expected, gradient of displacement from 
top to bottom shows that the system is under vertical 
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Fig. 9 Displacement field along the loading axis, (a) The overall 
vertical displacement field is represented very well by an expo- 
nential equation. For visual convenience every 30 data points is 
shown, (b) Longitudinal displacement for the 1st stage and the 
overall displacement. Here we choose a bin size of 0.4mm for av- 
eraging, (c) Average radial displacement of grains measured from 
the distance between the centre of packing its outer boundary. 



compression. The decrease of the slope with depth also 
indicates that part of this stress is screened by the time 
it reaches the bottom. A typical explanation of this is 
the Janssen effect [34], suggesting that friction on the 
lateral walls might be significant. To confirm this, we 
studied the vertical displacement profile in the radial 
direction (figure |9jb)) which confirms that grain move- 
ments closer to the outer region of cylinder are smaller 
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Fig. 10 Response of a single rubber ball to loading. 

than that of the central region. Therefore a shear com- 
ponent in the (z, r) direction is expected. Another fea- 
ture of the displacement field is the weak but noticeable 
radial displacement of the grain towards the boundary 
(see figure [9^ c)). This radial displacement can be clearly 
seen in figure [8ja). 

2.2.2 Contact mechanics 

In addition to the direct measurement of individual 
grain displacements, a number of other mechanically 
relevant features can be extracted from the tomogram. 
Upon compression, we observed that the distance be- 
tween grains varies, as well as shape and contact areas 
between grains. These quantities can a priori be used 
to quantify the force transmitted between the grains. In 
this section, we describe the implementation of a con- 
tact model that is suitable for stress calculations from 
3D images. 

Contact mechanics between solid objects have been 
studied for over a century. The seminal works of Hertz 
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Fig. 11 Solid lines: the mean vertical stress as a function of 
contact area offset. Dotted lines: the applied load, 740 Pa for 
stage 1, and 5000 Pa for stage 2. 



[35] and Mindlin [36 about the contact mechanics of 
solid bodies have provided analytical solutions for ide- 
alised cases, summarised in [37]. The contact force be- 
tween two elastic spheres can be calculated from the 
nonlinear Hertz-Mindlin model. The expression of the 
normal force between two contacting elastic spherical 
grains m and n with uncompressed radii R m and R n , 
made of the same material, is given by: 

_ 2 4G pi/ 2 >3/2 m 

Jmn — o i Wn' V-V 



31 



where R is the geometric mean of R m and R ni R = 
2R m R n /(R m + Rn), £ mn is the normal overlap, or pen- 
etration length, as depicted in Fig. [TO^ i, G is the shear 
modulus and v the Poisson ratio of the grains, equal to 
0.5 for rubber [2]. 

We have tested the mechanical response of our grains 
by compressing a few individual grains between two 
steel plates and measuring the force required as a func- 



tion of the gap between the plates (Fig. 10). The re- 
sponse of the grains is consistent with the Hertz model, 
and we have extracted from these graphs the value of 
the grains shear modulus, 850/cPa, which is a reason- 
able value for a commercial rubber. 

The Hertz-Mindlin theory therefore provides us with 
a suitable model to calculate the contact force from 
the grain geometry, as long as we have a good esti- 
mate of the normal overlap £ for each contact. Two 
options can be considered at this stage, i) The nor- 
mal overlap can be estimated from the distance be- 
tween grain centres. If the grains are located at po- 
sitions r m and r n , assuming grains remain spherical, 
£mn = \[Rm + Rn — |l m — r n |]. The autocorrelation 
function shows that we expect this quantity to be of the 
order of a few percent of the grain diameter. Although 
the grain location is determined accurately, calculation 
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Fig. 12 Histograms of normal forces in the sample for all three 
loading stages. 



of £ can be fairly inaccurate due to the anisotropy of 
the imposed strain causing the grains to deviate sig- 
nificantly from a spherical shape into more ellipsoidal 
geometries. The centre-centre distance is therefore not 
a good approach to estimate the contact geometry, ii) 
Our image processing protocols deliver a sensitive mea- 
surement of the contact area (Fig. [3]), which provides 
us with a more reliable way to calculate the forces in- 
dependently of the centre-centre distances. The contact 
area (s mn ), between grains labelled by m and n, is re- 
lated with the normal overlap £ mn by: £ mn = |™r . This 
purely local measurement of the force can be applied to 
any grain geometry, as long as the local curvature of the 
grain near the contact point can be estimated as well. 

2.2.3 Force distribution and stress field 

As discussed above, the finite resolution of the CT causes 
the contact areas to be systematically larger than their 
real values so this measurement needs to be carefully 
calibrated. There is de facto a small distance 5, of the 
order of the voxel size, such that if two surfaces are sep- 
arated by less than (5, they will appear in contact after 
segmentation. This results in a systematic enlargement 
of the contact area, and even detection of false con- 
tacts if the separation between grains is less than S. In 
the case of spherical elastic grains, the contact area af- 
ter segmentation would correspond to the geometrical 
overlap of two grains with each radius being increased 
by 6/ 2. Since the contact area is linear with the con- 
tact overlap £, the effect that S has on the surface area 
is equivalent to a systematic offset s Q ff = ttR5. How- 
ever, S is a priori unknown and has to be determined 
by calibrating our measurements. 

The mean stress in the vertical direction at the vicin- 
ity of the upper interface can be estimated in two in- 
dependent ways for each loading step; either from the 



knowledge of the loading mass and geometry of the 
setup, or by using the contact forces to calculate the 
stress tensor in the bulk of the pile. We consider a sub- 
volume Sv. Each contact contained in this volume bears 
a force denoted f mn , measured from the contact area 
Smn between grains indexed by m and n. The centre 
to centre vector is noted l m n . The components of the 
stress tensor, indexed by i and j , are then obtained from 
the following expression: 



(2) 



where the sum is over all contacting pairs of grains in 
the volume Sv. 

In order to measure the offset S required to com- 
pensate for the segmentation error, we have calculated, 
using the upper third of the sample, the mean vertical 
stress cr zz we would obtain for a range of values of the 
surface offset s Q ff 



ttRS (Fig. 11). A suitable choice 



should provide the value of the normal stress consistent 
with the loading applied on the sample (i.e. where the 



solid line intersects the dotted line on Fig. 11). Based 
on this graph, we select s Q ff « 0.25 mm 2 , that nearly 
sets the normal stress to its expected value for the high 
load value, where the calibration is the most reliable 
due to the large values of the force and increased num- 
ber of contacts. This value corresponds, as expected, to 
the size of a single voxel, confirming the consistency of 
the method. For stage and 1, however, based on this 
single threshold calibration, the overall mean stress ap- 
pears to have been overestimated. 

Once this tool is calibrated, it can be used to probe 
a number of statistical quantities in the pile. Figure 12 
shows the distribution of normal forces in the pile for 
the three different loads studied here. The 3D bulk mea- 
surement of the contact normal forces exhibit a number 
of features that are characteristic of granular systems. 
In particular, a salient feature of force measurement 
presented in this study is that the distributions are pri- 
marily exponential for large forces. This is in agree- 
ment with earlier measurements in 2D piles [38,39] or 
at a 3D interface [21] and also the numerical studies 
of dynamics of granulated systems [5] . Another feature 
often reported in granular dynamics studies [3|[4Q] is 
the presence of a peak or a plateau at low forces. How- 
ever, this region of the distribution is also where the 
forces are the most sensitive to image resolution, seg- 
mentation and thresholding. We believe a much higher 
spatial resolution is required before conclusions can be 
made about the precise shape of the force distribution 
in the region of low forces, in particular since the mean 
stress could not be adjusted. 
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Fig. 13 Stress fields in a plane section of the sample for various 
loads. From top to bottom, i) cartoons showing the location of 
the section in the sample, along a diameter, and various stress 
components reported, ii) the stress map in the azimuthal direc- 
tion (Jqq, in the radial direction a rr , in the vertical direction a zz , 
and shear stress o rz - 



The stress field in the pile is represented by a se- 
ries of color maps in Fig. [13] These stress maps con- 
firm a number of expected results. In the absence of 
any loading (stage 0), we observe a slight increase of 
the stress from top to bottom, in agreement with the 
fact that only the weight of the grains themselves acts 
at this stage. The same pattern is observed in all di- 
rections. The vertical component of the stress tensor 
shows an increase with the loading, as well as the ra- 
dial and azimuthal components. The magnitude of the 
vertical component is about twice as large as the other 
two. Another finding of interest concerns the existence 
of the shear component (r,z). It is also notable that 
both vertical and radial components of the stress tensor 
have larger values near the boundary of the confining 



cylindrical cell at the 2nd stage of compression. Grains 
at the centre experience a reduced compression, due to 
their ability to move slightly towards the sides. Lateral 
grains are however under stronger radial and vertical 
displacement gradients. These maps are therefore in full 
agreement with the displacement field measured from 
the tracking of grains for high loadings (Fig[9|. 

Conclusion 

We have presented in this paper the first measurements 
of internal stresses in a dry granular system resolved at 
the single grain scale. We have used here a simple geom- 
etry that allows us to calibrate and validate our mea- 
surements. In particular, we are now able to quantify, 
for the same pile, a number of characteristic features. 
They include mainly the evolution of the coordination 
number, the 2-point correlation function, the internal 
displacement fields, force distributions and stress fields, 
all as a function of loading. Taken together, all these 
measurements will enable us to better unravel the mi- 
cromechanical behaviour of granular systems, in partic- 
ular in quasistatic regimes. It is worth pointing out as 
well that by measuring the force based on the contact 
geometry, we are able to deal with non symmetric and 
polydisperse grains. 

Although we have proven that these measurements 
are realistic and achievable, a number of improvements 
are still required in order to study more complex cases. 
We need, in particular, to increase both the spatial 
statistics (number of grains) and the resolution at the 
contact scale which will provide a better calibration 
of the force and therefore a more reliable force mea- 
surement in the pile. This will also improve our mea- 
surements of coordination number and other geomet- 
ric quantities. Increased image resolution will also al- 
low the use of stiffer grains which in turn will widen 
the range of suitable materials to use for the grains. 
These improvements are achievable in the near future 
using high resolution nano-focus CT combined with 
large panel detectors, providing at least a five fold im- 
provement in resolution and speed. The contact model 
also needs further refinements, so that tangential forces 
can be estimated from the orientation of the contact 
zone with respect to the centre-to-centre line. These 
technical advances, combined with the analysis tools 
presented here, will contribute to the understanding of 
a number of open problems. Internal force distribution 
can now be studied in the bulk, reaching the low force 
region of the distribution. The mechanics of the mate- 
rial and in particular its mechanical stability rely on the 
internal organisation of the grain contacts and also the 
force network. These are quantities that can be directly 
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analysed in 3D, not only statically, but also dynami- 
cally by monitoring how an applied load affects each 
individual grain contact. 

The kinematic study of the material can also be 
qualitatively and quantitatively extended. Not only the 
local strain can be measured, but the non-affine part 
of the displacement field is also directly accessible. The 
latter is increasingly thought to be related with the bulk 
material properties, in particular when the system is at 
the onset of rigidity [41 . How such ideas can apply to 
3D frictional granular systems remains to be investi- 
gated. In particular, answering such questions will re- 
quire an extension to our current method so that we 
can study grain rotation. This would allow us to test 
the importance of micropolar elasticity in the mechan- 
ics of granular systems. 
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